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HIGHLIGHTS 


•  We  propose  an  enhanced  single  particle  (ESP)  model  for  lithium-ion  batteries. 

•  The  ESP  model  accounts  for  the  solution  phase  limitation. 

•  The  ESP  model  is  available  for  the  high  rate  charge-discharge  analysis. 

•  We  report  the  two-way  electrochemical-thermal  coupled  simulation  method. 

•  Temperature  estimates  during  the  charge-discharge  cycle  are  calculated  accurately. 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  11  October  2013 
Received  in  revised  form 
27  November  2013 
Accepted  29  November  2013 
Available  online  10  December  2013 


Keywords: 
Lithium-ion  battery 
Thermal  behavior 
Simulation 
Battery  modeling 
Single  particle  model 
Two-way  coupling 


To  understand  the  thermal  behavior  of  lithium-ion  secondary  batteries,  distributed  information  related 
to  local  heat  generation  across  the  entire  electrode  plane,  which  is  caused  by  the  electrochemical  re¬ 
action  that  results  from  lithium-ion  intercalation  or  deintercalation,  is  required.  To  accomplish  this,  we 
first  developed  an  enhanced  single  particle  (ESP)  model  for  lithium-ion  batteries  that  provides  a  cost 
effective,  timely,  and  accurate  method  for  estimating  the  local  heat  generation  rates  without  excessive 
computation  costs.  This  model  accounts  for  all  the  physical  processes,  including  the  solution  phase 
limitation.  Next,  a  two-way  electrochemical-thermal  coupled  simulation  method  was  established.  In  this 
method,  the  three  dimensional  (3D)  thermal  solver  is  coupled  with  the  quasi-3D  porous  electrode  solver 
that  is  applied  to  the  unrolled  plane  of  spirally  wound  electrodes,  which  allows  both  thermal  and 
electrochemical  behaviors  to  be  reproduced  simultaneously  at  every  computational  time-step.  The  quasi- 
3D  porous  electrode  solver  implements  the  ESP  model. 

This  two-way  coupled  simulation  method  was  applied  to  a  thermal  behavior  analysis  of  18650-type 
lithium-ion  cells  where  it  was  found  that  temperature  estimates  of  the  electrode  interior  and  on  the 
cell  can  wall  obtained  via  the  ESP  model  were  in  good  agreement  with  actual  experimental 
measurements. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Lithium-ion  batteries  are  attractive  power  sources  for  hybrid 
electric  vehicles  (HEVs)  and  electric  vehicles  (EVs)  because  they 
offer  high  cell  voltages  and  energy  densities.  As  a  result,  significant 
amounts  of  research  have  been  focused  on  improving  the  power 
density,  energy  density,  and  life  cycle  performance  of  such  batte¬ 
ries.  However,  because  one  of  the  primary  issues  in  designing 
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lithium-ion  batteries  is  ensuring  safety,  ceaseless  efforts  have  been 
expended  towards  overcoming  thermal  stability  problems. 

When  evaluating  thermal  stability  in  the  battery  design  phase, 
numerical  simulation  techniques  provide  useful  information  that  is 
difficult  or  impossible  to  obtain  experimentally.  As  a  result,  there 
have  been  numerous  previous  reports  on  thermal  modeling  of 
lithium-ion  batteries  over  a  wide  range  of  exposed  conditions. 
Many  such  reports  go  beyond  normal  use  [1—11  ,  with  some 
focusing  specifically  on  thermal  abuse  [12,13]. 

For  example,  Pals  and  Newman  performed  one-dimensional 
(ID)  thermal  modeling  experiments  to  calculate  temperature  pro¬ 
files  in  cell  stacks  [5,6  .  This  work  was  based  on  the  ID  macroscopic 
model  developed  by  Doyle  and  Newman  14]  with  the  addition  of  a 
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lumped  heat  generation  term  presented  by  Bernardi  et  al.  15]  Chen 
and  Evans  [7-9]  presented  a  multi-dimensional  thermal  model  for 
lithium-ion  batteries  that  focused  on  heat  transport  inside  the  cell 
stack  without  considering  the  electrochemistry  of  the  cell.  In  their 
work,  the  heat  generation  rate  was  estimated  using  experimental 
discharged  curves  based  on  the  lumped  heat  generation  formulas 
given  by  Bernardi  et  al.  Kim  et  al.  10,11  applied  a  two-dimensional 
(2D)  model  to  parallel  plate  battery  electrodes  in  order  to  simulate 
not  only  the  potential  and  current  density  distribution,  but  also  the 
temperature  distribution.  However,  their  modeling  approach 
focused  on  ensuring  current  continuity  on  the  electrodes  and 
lumped  heat  generation  terms  in  a  fashion  similar  to  the  above- 
mentioned  models. 

In  this  study,  we  focus  on  the  thermal  behaviors  of  lithium-ion 
batteries  under  normal  operating  conditions.  The  objective  of  our 
work  is  the  development  of  a  multi-dimensional  simulation 
method  capable  of  evaluating  both  electrochemical  and  thermal 
behaviors  inside  lithium-ion  batteries  as  accurately  as  possible, 
without  the  use  of  a  lumped  heat  generation  formula. 

In  order  to  achieve  this,  it  was  first  necessary  to  develop  a  new 
variation  of  the  commonly  used  single  particle  (SP)  model  for 
lithium-ion  batteries.  The  previous  simplified  battery  model  [16- 
21]  neglects  the  solution  phase  limitation.  In  contrast,  the  new 
and  simplified  battery  model  presented  in  this  paper  is  capable  of 
estimating  heat  generation  rates  using  a  detailed  theoretical  for¬ 
mula  [1-3],  since  it  accounts  for  all  the  physical  processes, 
including  the  solution  phase  limitation. 

Subsequently,  a  quasi-3D  porous  electrode  solver  was  devel¬ 
oped  implementing  the  above  new  simplified  battery  model.  This 
solver  is  basically  applied  to  the  2D  unrolled  plane  of  spirally 
wound  electrodes  and  provides  distributed  information  on  the 
potentials  and  lithium-ion  concentrations,  not  only  across  the 
entire  2D  electrode  plane,  but  also  across  the  thickness  of  the 
electrode  itself,  which  can  be  obtained  with  the  assistance  of  the 
new  simplified  battery  model.  As  a  result,  the  local  heat  generation 
rate  resulting  from  the  electrochemical  reactions  that  occur  during 
the  charge-discharge  process  can  be  calculated  in  accordance  with 
the  detailed  theoretical  formula,  and  its  2D  distribution  data  can  be 
obtained  across  the  entire  electrode  plane. 

Finally,  a  two-way  electrochemical-thermal  coupled  simulation 
method  has  also  been  established.  Rising  temperatures  in  a 
lithium-ion  battery  affect  its  performance  because  both  the  trans¬ 
port  and  physical  properties  strongly  depend  on  temperature.  In 
this  method,  the  quasi-3D  porous  electrode  solver  transfers  the 
information  of  the  planer  distribution  of  the  heat  generation  rate  to 
the  3D  thermal  solver.  Following  that,  the  3D  thermal  solver  sim¬ 
ulates  the  temperature  distribution  of  the  spirally  wound  elec¬ 
trodes  and  returns  the  temperature  information  to  the  quasi-3D 
porous  electrode  solver.  By  executing  the  above-mentioned  data 
exchange  process,  the  temperature  dependency  of  the  transport 
and  physical  properties  can  be  considered  at  every  computational 
time-step. 

This  paper  describes  the  abovementioned  new  charge- 
discharge  SP  model  for  lithium-ion  batteries  and  the  two-way 
electrochemical-thermal  coupled  simulation  method  in  detail.  In 
Sections  2  and  3,  model  development  and  validation  of  the  new 
simplified  battery  model  are  shown  in  comparison  with  the  con¬ 
ventional  SP  model.  In  Section  4,  the  two-way  electrochemical- 
thermal  coupled  simulation  method  is  described.  In  Section  6, 
simulated  results  of  the  electrochemical  and  thermal  behaviors  of  a 
18650-type  lithium-ion  cell  during  the  charge-discharge  cycle  are 
expressed  and  discussed  in  detail  with  the  aid  of  temperature 
measurements  obtained  from  inside  a  real  cell,  which  is  shown  in 
Section  5.  Finally,  conclusions  of  this  paper  are  presented  in 
Section  7. 


2.  Model  development 

In  this  section,  the  modeling  approach  for  our  new  charge- 
discharge  SP  model  for  lithium-ion  batteries  is  described  in  detail. 

2.1.  Conventional  lithium-ion  battery  models 

The  literature  on  modeling  approaches  for  lithium-ion  batteries 
is  quite  extensive.  However,  conventional  mathematical  models 
can  be  classified  roughly  into  two  categories. 

The  first  modeling  approach  was  proposed  by  Doyle  and  New¬ 
man  14].  A  schematic  of  this  model,  which  consists  of  two  com¬ 
posite  electrodes  and  a  separator,  is  shown  in  Fig.  1  and  is  referred 
to  as  the  Newman  model  hereafter.  A  mathematical  representation 
includes  equations  that  describe:  (1)  mass  transport  of  lithium  in 
the  solid  phases,  (2)  mass  transport  of  lithium-ions  in  the  solution 
phase,  (3)  charge  transport  in  the  solid  phases,  and  (4)  charge 
transport  in  the  solution  phase.  The  Newman  model  is  treated  as  a 
ID  macroscopic  model  across  the  thickness  of  the  electrode  in 
the  local  point  on  the  electrode  plane.  To  put  it  more  precisely, 
the  concentration  of  lithium  within  the  solid  phase  is  solved 
rigorously,  using  the  extra  pseudo  second  dimension  along  the 
radius  of  the  particle.  From  the  viewpoint  of  application  to  the  3D 
electrochemical-thermal  coupled  simulation  for  lithium-ion  bat¬ 
teries,  the  advantage  of  the  Newman  model  is  its  ability  to  accu¬ 
rately  estimate  heat  generation  rates,  whereas  its  demerit  is  its  high 
computational  costs.  The  local  heat  generation  rate  is  calculated 
from  the  detailed  theoretical  formula  shown  below  [1-3]. 

q  =  a  effV0s-V0s+  ^xeffV0e  •  V0e  +  In  ce  •  +  as,  jin,  j  (fis 

9  Uj 

-fie  ~  Uj)  +  as,jln,jT-QY 

(1) 

Eq.  (1 )  requires  the  local  point  values  and  gradients  of  potentials 
and  concentrations  in  the  solid  and  electrolyte  phase.  Since  the 
Newman  model  can  be  used  to  estimate  these  values  at  the  each 
local  point,  the  heat  generation  rate  can  be  accurately  estimated  via 
Eq.  (1 ).  On  the  other  hand,  the  Newman  model  requires  the  use  of  a 
pm-order  mesh  size  to  discretize  a  calculation  domain  across  the 
thickness  of  each  electrode  in  ID.  As  a  result,  if  the  entire  lithium- 
ion  cell  geometry  of  a  18650-type  cell  is  discretized  with  the  pm- 
order  size,  the  number  of  computational  meshes  will  be  far  in 
excess  of  ten  million,  which  makes  this  methodology  expensive  and 
time  consuming. 


Anode  Cathode 


Fig.  1.  Schematic  of  the  lithium-ion  battery  model  proposed  by  Doyle  and  Newman 

[12]. 
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The  second  modeling  approach  is  the  SP  model  [16-21  .  A 
schematic  of  this  model  is  shown  in  Fig.  2.  In  this  model,  each 
electrode  is  represented  by  a  single  spherical  particle.  When 
considering  the  3D  electrochemical-thermal  coupled  simulation  of 
lithium-ion  batteries,  it  is  important  to  reflect  on  the  advantages 
and  drawbacks  of  the  SP  model.  The  advantages  of  the  SP  model 
include  its  simplicity  and  reduced  computation  time  requirements. 
This  model  is  orders  of  magnitude  faster  than  the  Newman  model. 
However,  a  significant  drawback  of  the  SP  model  is  that  it  assumes 
that  the  electrolyte  phase  diffusion  limitations  are  ignored.  Hence, 
potentials  in  the  electrolyte  phase  are  set  to  zero,  and  concentra¬ 
tions  in  the  electrolyte  phase  are  set  to  a  constant  value.  This  ren¬ 
ders  estimations  of  the  heat  generation  rate  via  Eq.  (1)  inadequate. 
Furthermore,  even  if  the  heat  generation  rate  is  estimated  using  the 
lumped  formula  presented  by  Bernardi  et  al.  15],  the  resulting 
estimated  values  seem  to  have  less  accuracy. 


Anode  Cathode 


Fig.  3.  Schematic  of  the  ESP  model.  In  this  model,  each  of  negative  and  positive 
electrodes  is  represented  by  a  single  spherical  particle  in  the  electrolyte  phase. 


2.2.  Description  of  new  SP  model 

In  the  above  discussion,  the  advantages  and  drawbacks  of  con¬ 
ventional  lithium-ion  battery  models  were  classified.  In  the  present 
paper,  from  the  viewpoint  of  application  to  the  3D  electrochemical- 
thermal  coupled  simulation,  we  propose  a  new  lumped  model  for 
lithium-ion  batteries. 

This  model  combines  advantages  of  both  the  Newman  and  SP 
models.  However,  to  be  successful,  the  new  lumped  model  must 
provide  the  following  advantages: 

(a)  Calculation  costs  equivalent  to  that  of  the  SP  model 

(b)  Heat  release  rate  estimation  accuracy  at  same  level  as  the 
Newman  model 


model.  Mass  transport  and  charge  transport  between  both  elec¬ 
trodes  in  the  electrolyte  phase,  which  satisfies  the  mass  and  charge 
conservation,  are  estimated  using  these  boundary  values. 


2.3.  Model  equations 

23 A.  Calculation  procedure  of  concentrations 

Fig.  4  describes  the  procedure  for  calculating  ESP  model  con¬ 
centrations  in  detail. 

Fields  second  law  usually  describes  the  lithium  transport  in  the 
solid  phase,  with  the  following  boundary  conditions  set  for  a 
spherical  particle: 


The  development  concept  of  a  new  lumped  model  is  based  on 
the  SP  model  with  specific  consideration  paid  to  the  electrolyte 
phase  diffusion  limitations.  Fig.  3  shows  a  schematic  of  a  new 
lumped  model,  which  is  referred  to  hereafter  as  the  enhanced 
single  particle  (ESP)  model.  In  the  ESP  model,  each  electrode  is 
represented  by  a  single  spherical  particle  along  with  its  electrolyte 
phase.  The  potential  and  lithium  concentration  in  the  solid  phase 
are  calculated  in  the  same  way  as  with  the  conventional  SP  model. 
Additionally,  the  potential  and  lithium-ion  concentration  in  the 
electrolyte  phase  are  set  at  the  representative  position  in  the  each 
electrode. 

The  most  important  advancement  over  the  SP  model  at  this 
point  is  that  these  physical  quantities  are  approximated  by  para¬ 
bolic  profiles  within  each  electrode.  As  described  in  Eq.  (1 ),  accurate 
heat  release  rate  estimates  require  the  point  values  and  gradients  of 
these  physical  quantities.  These  point  values  and  gradients  are 
calculated  from  approximated  parabolic  profiles,  which  makes  it 
possible  for  the  ESP  model  to  calculate  the  heat  release  rate. 
Furthermore,  the  potentials  and  lithium-ion  concentrations  at  in¬ 
terfaces  between  the  negative  and  the  separator,  or  between  the 
positive  and  the  separator,  are  implicitly  considered  in  the  ESP 
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Fig.  2.  Schematic  of  the  single  particle  (SP)  model.  This  model  assumes  that  the 
electrolyte  phase  diffusion  limitations  are  ignored. 


dCs  „  d2Cs  2dCs 
at  s  ar2  r  dr 


(2) 


0  at  r  =  0 


(3) 


in 

F 


at  r  =  rs 


(4) 


In  the  conventional  SP  model,  the  diffusion  length  method  [22] 
is  introduced  to  simplify  the  diffusion  equation.  By  assuming  a 
parabolic  concentration  profile  in  the  diffusion  layer  and  using  the 
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volume  average  technique,  the  solutions  for  Eqs.  (2)— (4)  can  be 
approximated  with  a  set  of  differential  and  algebraic  equations: 


d csj  1 5D5  j 

dt  r2  ■ 

S,J 


0 


(5) 


asJF 


(6) 


where  lsej  is  the  diffusion  length  and  takes  the  value  of  rSJ/5  for 
spherical  particles.  In  view  of  the  need  to  reduce  computational 
costs,  the  diffusion  length  method  is  implemented  to  the  ESP  model 
in  the  same  way. 

The  diffusion  limitations  in  the  electrolyte  phase  are  also 
considered  in  the  ESP  model.  The  representative  lithium-ion  con¬ 
centrations  of  electrolyte  in  each  electrode  region  are  defined  by 
using  the  volume  average  technique  as  follows: 


1  Ln 

Ce,  n  =  j—  J  Q,  ndZn 


(7) 


Q,  p 


1 

Lp 


Q,  pdZp 


(8) 


In  the  ESP  model,  a  parabolic  concentration  profile  is  assumed  in 
each  region.  Therefore,  the  positions  corresponding  to  the  repre¬ 
sentative  concentrations  are 


* 

7 

^n.  e 


(9) 


* 

ZP-  e 


(10) 


where  the  extra  pseudo  coordinates  zn  and  zv  are  defined  as  shown 

in  Fig.  4. 

In  addition,  as  Fig.  4  also  shows,  the  concentrations  on  the 
negative/separator  interface  and  the  separator/positive  interface 
are  implicitly  defined  as  \  and  Qf  2,  respectively.  The  interfacial 
balances  of  lithium-ion  flux  yield  are 

neff  ^if ,  1  Ce,  n  neff  2  ^  fl  1  ^ 

ue,  n  X  ~  ue,  sep  r  v 1  1 ; 


neff  Cif ,  2  “  C(M  _  neff  ce ,  P  ~  cif ,  2 

ue.  sep  r  ~  ue,p  * 

Csep  Op 


(12) 


where  5n  and  dp  represent  the  diffusion  length  in  each  region  which 
are  expressed  as  [22] 


<5„  =  y  and  8P  =  ^ 


(13) 


Solving  Eqs.  (11)  and  (12)  gives  the  interfacial  concentrations  as 


Ci/,  1  — 


Oiyi  H“  CO 


ocn  +  co  +  y  p 


■  Q,  n  + 


vP 


0Cn  +  CO  +  Jp 


’ce,  p 


(14) 


Cif.  2  - 


(Xn 


OCn  +  co  +  Jp 


■  Q,  n  + 


w  +  Tp 


+  0)  +  Yp 


p 


(15) 
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Fig.  5.  Calculation  procedure  of  the  potential  field  of  the  ESP  model. 
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The  equations  for  calculating  the  time  evolution  of  the  volume 
averaged  lithium-ion  concentration  in  the  electrolyte  phase  are 
derived  as  shown  below.  As  for  the  negative,  conservation  of 
lithium-ion  [14]  in  the  electrolyte  yields 


0Q,  n  _  neff  ce,  n  1  -Li 

fe-  n~dT  ~  Ue -  n_az2~~ +  ~T~Jn 


(17) 


As  shown  in  Fig.  4,  the  electrolyte  phase  concentration  is  rep¬ 
resented  by  a  parabolic  profile: 


Ce,  n  —  Cl’Z^  b-Zn  C  (18) 

where  the  three  coefficients,  a,  b  and  c  can  be  determined  by  the 
following  boundary  conditions  and  Eq.  (7). 

"  =  0  at  zn  =0  (19) 

dZn 


Ce,  n  —  Cif  -j  at  Zn  —  Ln  (20) 

Application  of  volume  averaging  to  Eq.  (17)  and  substitution  of 
Eq.  (18)  results  in 


£e,  n 


dCe,  n 

dt 


D 


eff 
e,  n 


(2a)  + 


1  -t° 


7  Li 
Jn 


Substituting  the  coefficient  a  into  Eq.  (21)  yields 


(21) 


dc, 


Ce,  rr 


e,  n 


dt 


—  ’  Ce,  n  Bn  -  Ce?  p  + 


1  -t° 


Li 


Jn 


(22) 


where  the  coefficients  An  and  Bn  are  given  by 


J_  Jp'OCn 
Bn  CXn  +  CO  +  J p 


(23) 


1  Tp-^n 
£71  U!n  +  CO  +  Jp 


(24) 


As  for  the  positive,  the  same  procedure  mentioned  above  also 
yields  the  equation  for  calculating  the  time  evolution  of  the  volume 
averaged  lithium-ion  concentration  in  the  electrolyte  phase: 
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d  Ce  n  A  n  _  1  t,  .i  ; 

£e,  P  fa  =  Ap'ce,n  +  Bp-  Ce,  p  H  - — Jp 


where  the  coefficients  Ap  and  Bp  are  given  by 

A  =  1.  JP'an 
p  Lp  otn  +  (o  +  yp 


(25) 


«app 

Bn 


—  as,  rd 0,  n 


f  &a,  nB 

6XP(  RT  V'[ 


exp 


«c,  nB 

RT 


Vn 


(35) 


where  the  surface  overpotential,  rjnt  is  defined  to  be 


(26)  Vn  —  4s,  n  4e,  n  Un  ^n,n‘Byn  (56) 

Conservation  of  charge  in  the  electrolyte  phase  results  in  14] 


Bp  =  -f - —  (27) 

Bp  Oin  +  (0  +  Yp 

Currently,  in  the  ESP  model,  Eqs.  (5),  (6),  (22)  and  (25)  are 
specified  in  order  to  solve  the  mass  transport  of  the  system. 

2.3.2.  Calculation  procedure  of  potentials 

With  the  assistance  of  Fig.  5,  the  calculation  procedure  of  po¬ 
tentials  of  the  ESP  model  is  described  in  detail.  In  analogy  with  the 
calculation  procedure  of  concentrations,  the  potentials  on  the 
negative/separator  interface  and  the  separator/positive  interface 
are  implicitly  defined  as  fry  i  and  fry  2,  respectively. 

The  boundary  conditions  for  the  solid  phase  potential  0s>n  in  the 
negative  electrode  are  described  below.  On  the  copper  current 
collector/composite  negative  electrode  interface,  fr<n  is  arbitrarily 
set  to  zero,  while  on  the  same  interface,  the  solid  phase  current 
density  is  equated  to  the  applied  current  density,  namely 


4>s,  n  =  0  at  zn  =  0 

(28) 

eff9<^  n  _  ,  ,t  z  ~  0 

°n  qz  —  Upp  ai  Zn  —  U 

(29) 

On  the  negative/separator  interface,  the  charge  flux  in  the  solid 
phase  is  equated  to  zero. 

dtl’n  =0  at  Zn  =ln 

9 Zn 

(30) 

Applying  a  parabolic  profile  for  0s,m  satisfaction  of  the  above 
boundary  conditions  results  in 

A  1  ^pp  ^2  Japp  ^ 

n  ~  2 Bn  (7eff  (7eff  " 

On  0  n 

(31) 

The  representative  potential  in  the  negative 
defined  by  using  the  volume  average  technique. 

solid  phase  is 

1  Ln 

4s,  n  =  J  n^Zn 

0 

(32) 

V-  (k fvtej)  +  v-  feDffjV  In  Cej)  +jf  =  0  (37) 

From  Eq.  (37),  the  charge  flux  on  the  negative/separator  inter¬ 
face  is  equated  to  the  applied  current  density: 


—K 


:9</> 


9c, 


e,  n 


eff  Ye,  n 


n 


dZn 


-  K 


Zn — tn 


eff 
D,  n 


dzn 


Zn — tn 


Cif ,  1 


fapp 


(38) 


Since  the  second  term  of  the  left-hand  side  is  estimated  using 
Eqs.  (14)  and  (18),  the  potential  gradient  on  the  negative/separator 
interface  (the  first  term)  is  obtained  from  Eq.  (38).  Next,  a  parabolic 
profile  is  assumed  for  the  electrolyte  potential  across  the  thickness 
of  the  negative  electrode.  This  parabolic  profile  must  meet  the 
following  three  conditions: 


®4e,  n 

dZn 


at  zn  =  0 


(39) 


®4e,  n 

dZn 


sif,  1  at  zn 


Bn 


(40) 


4e,  n  —  ^  J  4e,n^zn  (41) 

0 

where  Sy  i  is  the  electrolyte  potential  gradient  on  the  negative/ 
separator  interface,  which  is  estimated  from  Eq.  (38).  Using  Eqs. 
(39)— (41),  the  potential  profile  in  the  negative  electrolyte  region  is 
found  to  be  determined  by 


4e,  n 


Sif,  1 
2Bn 


-Zn  +  4e,  n 


(42) 


According  to  Eq.  (42),  the  interfacial  potential  fry  \  is  readily 
found  to  be  expressed  by 


4* if,  1  —  — Bn  +  4e,  n  (43) 

Taking  into  account  the  continuity  of  the  charge  flux  in  the 
electrolyte  phase,  the  charge  transport  in  the  separator  leads  to 


Substituting  Eq.  (31)  into  Eq.  (32)  gives 


4s,  n 


1  Bn 

3  oW 


(33) 


This  is  equal  to  the  point  value  that  is  calculated  in  the  profile, 
Eq.  (31),  at 


* 

zn,  s 


34 


4 if,  2  —  4 if,  1  -  ^eff  ^app-^sep  +  ^sep  (44) 

where  the  interfacial  concentrations,  qt  i  and  C\r  2,  can  be  estimated 
by  Eqs.  (14)  and  (15). 

A  similar  calculation  procedure,  which  is  mentioned  above, 
derives  the  parabolic  profile  for  the  positive  electrolyte  potential. 
The  charge  flux  on  the  positive/separator  interface  is  equated  to  the 
applied  current  density: 


The  calculation  procedure  of  the  other  potentials  (0e  n,  fry  7,  fry  2, 
4e,  p>4s,  p)  are  described  sequentially  below. 

Because  Eq.  (33)  gives  the  solid  phase  potential,  </>s  n,  the 
representative  potential  in  the  negative  electrolyte  region,  </>e  n,  is 
determined  by  the  following  Butler-Volmer  equation,  namely 
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:d(f> 
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eff  Ye,  p 
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zp= 0 


eff 
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dz 
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cif.2 
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app 


(45) 


The  second  term  on  the  left  is  estimated  by  using  Eq.  (15)  and 
the  parabolic  profile  of  ce>  p.  This  allows  the  potential  gradient  on 
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the  separator/positive  interface  in  the  first  term  to  be  estimated.  A 
parabolic  profile  is  assumed  for  the  electrolyte  potential  across  the 
thickness  of  the  positive  electrode.  This  parabolic  profile  must 
satisfy  the  following  conditions: 


0e,  p  ~~  fiif,  2  at  Zp  —  0 

(46) 

90e  p  n 

dZpP=siL2  at  zp  =  0 

(47) 

90e,  p  n  T 

_0  at  zp  -Lp 

(48) 

2.3.3.  Summary  of  model  equations 

Model  equations  of  the  ESP  model  are  now  specified  completely. 
The  final  equations  are  summarized  in  fable  1.  Time  evolutions  of 
unknowns  are  calculated  by  using  the  set  of  equations  for  volume 
averaging  values. 

Here,  the  relationship  between  the  conventional  SP  model  and 
the  ESP  model  is  worth  a  passing  mention.  The  conventional  SP 
model  assumes  that  the  diffusion  limitations  in  the  electrolyte 
phase  are  negligible.  This  assumption  leads  to  omitting  all  equa¬ 
tions  in  the  electrolyte  phase  from  the  ESP  model  presented  in  this 
paper.  The  conservation  of  charge  in  the  negative  solid  phase  is 
given  by 


where  s,£  2  is  the  electrolyte  potential  gradient  on  the  separator / 
positive  interface,  which  is  estimated  from  Eq.  (45).  Using  Eqs. 
(46)— (48),  the  potential  profile  in  the  positive  electrolyte  phase  is 
found  to  be  expressed  by 


9  0s,  n 

dz2 

U4n 


0 


Applying  Eq.  (31)  to  Eq.  (58)  yields 


(58) 


0e,  p  —  — ^jF — zp  +  sif,  2  %ZV  +  0(f,  2  (49) 

The  representative  potential  in  the  positive  electrolyte  is 
defined  to  be 


=  0  (59) 

Ln 

The  above  equation  is  also  used  in  the  conventional  SP  model. 
Note  that  the  ESP  model  that  neglects  the  electrolyte  phase  diffu¬ 
sion  limitations  reduces  to  the  conventional  SP  model. 


Lp 

0e,  p  =  J  0e,  pd-Zp  =  — ^  Lp  +  (fry  2  (50) 

P  0 

This  is  equal  to  the  point  value  that  is  calculated  in  the  profile  of 

Eq.  (49)  at  z*  e  (Eq.  (10)). 

Finally,  let  us  turn  to  the  positive  solid  phase  potential.  As  Eq. 
(50)  gives  the  electrolyte  phase  potential,  the  representative  po¬ 
tential  in  the  positive  solid  phase,  </>s  p,  is  estimated  by  the 
following  Butler- Volmer  equation,  namely 


as,  pio,  p 


exp 


( ^a,  pF 

v  RT 


exp 


OLc,  pF 


RT 


where  the  surface  overpotential,  rjp,  is  defined  to  be 


(51) 


Vp  —  0s,  p  0e,  p  Lfp  in,p-Rf,p  (52) 

Just  as  in  the  case  of  the  solid  phase  potential  in  the  negative 
electrode,  0s,m  a  parabolic  profile  is  assumed  and  must  meet  the 
following  conditions  to  determine  the  three  coefficients. 
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Satisfaction  of  above  conditions  results  in 
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2.3.4.  Estimation  of  heat  generation  rate 

The  heat  generation  rate  during  the  charge-discharge  process, 
which  is  Eq.  (1)  as  revised  for  the  ESP  model,  is  expressed  by 


V  {o fv<l>s,j-v<l>s,j)  +  E 

j  =  n,  p  j  =  n,  sep,  p 

+  KDJ^  1°  ce,  j '  V0e,  j)  +  ^2  a s>  j*n-  j  j  ~  j  ~  Uj) 


J  =  n,p 


~  dUj 

+  ^2  asJln:jT 

j  =  n,  p 


d  T 


(60) 


The  first  term  results  from  the  ohmic  heat  in  the  solid  phase, 
while  the  second  term  arises  from  the  ohmic  heat  in  the  electrolyte 
phase.  The  summation  of  the  last  two  terms  accounts  for  irre¬ 
versible  and  reversible  heats  associated  with  charge  transfer  at  the 
solid/electrolyte  interfaces.  The  third  term  represents  the  potential 
deviation  from  the  equilibrium  potential  (irreversible  heat).  The 
fourth  term  arises  from  the  entropic  effect  (reversible  heat). 

From  Eq.  (60),  it  is  clear  that  the  estimation  of  ohmic  heats  in  the 
solid  and  electrolyte  phases  needs  the  information  about  the  gra¬ 
dients  of  the  unknowns.  The  ESP  model  includes  the  profile  infor¬ 
mation  of  concentrations  and  potentials  across  the  thickness  of  the 
electrode  as  standard  features,  as  shown  in  Table  1.  The  gradients  of 
concentrations  and  potentials  in  the  solid  phase  are  estimated  in 
each  parabolic  profile  at  the  positions  of  z*  s  (Eq.  (34),  for  the 
negative)  and  z*  s  (Eq.  (57),  for  the  positive).  Additionally,  the 
gradients  in  the  electrolyte  phase  are  estimated  in  each  parabolic 
profile  at  the  positions  of  z*  e  (Eq.  (9),  for  the  negative)  and  z*  e  (Eq. 
(10),  for  the  positive).  As  for  estimation  of  irreversible  and  revers¬ 
ible  heats,  volume  averaging  values  are  used. 


3.  Model  validation 


The  representative  value,  </>sp,  is  obtained  in  the  profile,  Eq.  ,  3  7  Comparison  with  Newman  model  accuracy 

at 
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First  we  examined  the  cell  potential  during  galvanostatic 
(57)  discharge  at  various  current  densities  for  the  cell  described  in 
Ref.  23.  This  cell  consists  of  a  carbon  negative  electrode  (LixCe),  a 
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Table  1 

Summary  of  model  equations. 


Variable 
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in  solid  phase 
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j  =  n,  p 
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manganese  oxide  positive  electrode  (LiyMn204)  and  a  plasticized 
electrolyte.  The  solution  used  in  the  cell  is  2  M  LiPF6  in  a  1:2  ratio 
mixture  of  EC/DMC  in  the  plasticized  polymer  matrix.  The  1C  rate  of 
discharge  is  1.75  mA  cm-2.  For  the  values  of  the  parameters  used  in 
the  simulations,  see  Ref.  [23  . 

Fig.  6  shows  a  comparison  of  simulated  discharge  curves  be¬ 
tween  the  Newman  and  conventional  SP  models.  In  Fig.  6,  the  cell 
potential  is  given  as  a  function  of  the  capacity  (mAh  cm-2)  for 
discharges  at  various  rates.  As  pointed  out  in  Ref.  [17  ,  the  con¬ 
ventional  SP  model  is  as  valid  as  the  Newman  model  for  discharge 
rates  up  to  1C.  However,  the  slight  discrepancy  from  the  Newman 
model  observed  even  at  the  0.5C  rate  should  not  be  overlooked.  As 
stated  in  Section  2.1,  the  conventional  SP  model  assumes  the  con¬ 
centration  and  potential  in  the  electrolyte  phase  to  be  constant 
throughout  the  discharge.  Even  at  low  discharge  rates,  the  con¬ 
centration  change  in  the  electrolyte  phase  arises  moderately.  The 
conventional  SP  model  completely  ignores  this  change,  while  the 
Newman  model  follows  it  via  electrolyte  phase  charge/species 
conservation  equations.  This  is  reflected  in  the  above  slight 
discrepancy.  As  observed  from  the  figure,  the  difference  between 
the  two  models  becomes  increasingly  large  at  rates  of  discharge 
beyond  1C.  This  suggests  that  the  concentration  and  potential 
gradients  in  the  electrolyte  phase  limit  calculation  accuracy  under 
these  conditions. 

On  the  other  hand,  validation  of  the  ESP  model  presented  in  this 
paper  is  clearly  displayed  in  Fig.  7.  It  is  important  to  note  that  a 
significantly  good  degree  of  agreement  is  observed  between  the 


discharge  curves  of  both  the  Newman  and  ESP  models  for  rates  up 
to  2C.  The  ESP  model  is  capable  of  capturing  the  time  variations  of 
concentration  and  potential  in  the  solution  phase  during  the 
discharge  process  through  the  use  of  the  equations  shown  in 
Table  1.  This  advantage  over  the  conventional  SP  model  leads  to  a 
good  degree  of  agreement.  At  the  3C  or  4C  rate,  the  discharge  curve 
of  the  ESP  model  is  calculated  below  the  result  of  Newman  model  in 
the  early  discharge  stages,  in  contrast  with  the  results  of  the  con¬ 
ventional  SP  model  shown  in  Fig.  6.  As  a  result,  the  following 
additional  efforts  have  been  made  to  reduce  this  ESP  model  error. 

As  with  the  conventional  SP  model,  the  ESP  model  employs  the 
diffusion  length  method  inside  a  spherical  particle  with  a  second- 
degree  polynomial  (Eqs.  (5)  and  (6)).  The  method,  in  itself,  should 
be  valid  only  after  the  diffusion  layer  builds  up  to  its  steady  state. 
Therefore,  under  very  high  rates  of  discharge,  the  approximation 
that  fixes  the  diffusion  length  at  a  constant  steady  state  value  from 
the  beginning  of  the  discharge  is  inadequate  [20,24].  In  view  of  the 
method  shortcoming  mentioned  above,  the  time-variant  diffusion 
length  was  proposed  by  Wang  and  Srinivasan  [24]: 


exp 


(61) 


where  l*e  .  takes  the  value  of  rSJ/5  for  spherical  particles,  and  \sj  is 
the  tuning  parameter.  The  empirical  term  was  formed  based  on  the 
observation  that  the  surface  concentration  increases  exponentially 


Capacity  [mAh/cm2] 


Fig.  6.  Comparison  between  discharge  curves  at  various  rates  from  the  Newman 
model  and  those  predicted  by  the  conventional  SP  model. 


Capacity  [mAh/cm2] 


Fig.  7.  Comparison  between  discharge  curves  at  various  rates  from  the  Newman  model 
and  those  predicted  by  the  enhanced  SP  model  (ESP  model). 
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Fig.  8.  Comparison  between  discharge  curves  at  various  rates  from  the  Newman 
model  and  those  predicted  by  the  revised  ESP  model. 
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at  short  times.  Furthermore,  the  ESP  model  also  employs  the 
diffusion  length  method  concerning  the  electrolyte  phase  diffusion. 
The  diffusion  length  is  represented  by  Eq.  (13).  Based  on  the  similar 
discussion  as  in  the  case  of  the  solid  phase  diffusion,  the  time- 
variant  diffusion  length  is  introduced  to  the  ESP  model. 


(62) 


Here,  6*  is  given  by  Eq.  (13),  and  Xej  is  the  tuning  parameter. 

Fig.  8  shows  a  comparison  between  the  simulated  discharge 
curves  of  the  Newman  and  revised  ESP  models  using  Eqs.  (61)  and 
(62).  These  simulation  results  clearly  indicate  that  the  revised  ESP 
model  can  provide  results  that  are  as  accurate  as  the  Newman 
model,  even  at  high  discharge  rates. 


3.2.  Experimental  verification 

In  this  section,  experimental  verification  of  the  ESP  model  is 
described.  The  comparison  of  the  simulated  discharge  curves  for 
the  18650-type  lithium-ion  cell  is  shown  in  Fig.  9.  This  cell  was 


Distance  from  the  Cu  Current  Collector  [pm] 

Fig.  10.  Validation  of  the  predicted  profiles  in  the  electrolyte  phase  concentration  and 
potential  across  the  thickness  of  the  electrode  at  DOD  =  30%  during  the  5C  discharge 
rate  (in-house  manufactured  18650-type  lithium-ion  cell).  The  solid  lines  denote 
model  predictions  by  the  ESP  model  while  the  dashed  lines  denote  the  Newman  model 
predictions. 


manufactured  in-house  and  consists  of  a  carbon  negative  electrode, 
a  Ni-based  oxide  positive  electrode  and  a  separator.  The  capacity  of 
the  cell  is  approximately  0.5  Ah.  As  observed  at  a  2C  discharge  rate, 
both  the  conventional  SP  and  ESP  models  show  good  agreement 
with  the  experimental  discharge  curve.  However,  there  is  a  signif¬ 
icant  difference  between  simulated  results  under  very  high 
discharge  rates.  At  the  5C  and  10C  discharge  rates,  both  the  New¬ 
man  and  ESP  models  show  good  degrees  of  agreement  with  the 
experimental  discharge  curves.  More  specifically,  the  ESP  model  is 
able  to  provide  predictions  that  are  as  accurate  as  the  Newman 
model.  On  the  other  hand,  the  simulated  discharge  curves  obtained 
via  the  conventional  SP  model  show  a  significant  discrepancy 
beyond  the  experimental  results,  which  denote  the  same  tendency 
observed  in  Fig.  6. 

Incidentally,  the  ESP  model  provides  the  profile  information  for 
concentrations  and  potentials  across  the  thickness  of  the  electrode, 
as  shown  in  Table  1.  Fig.  10  shows  a  comparison  of  the  concentra¬ 
tion  and  potential  profiles  in  the  electrolyte  phase  across  the 
thickness  of  the  electrode  when  the  depth  of  discharge  (DOD) 
equals  30%  during  the  5C  discharge  rate.  This  figure  shows  that 
concentration  and  potential  profiles  across  the  thickness  can  be 
simulated  by  the  ESP  model  as  accurately  as  by  Newman  model. 

Finally,  let’s  summarize  the  advantages  of  the  ESP  model  pre¬ 
sented  in  this  paper. 

(a)  The  ESP  model  is  based  on  the  same  mass-point  system  used 
by  the  conventional  lumped  model.  Therefore,  this  model  is 
inexpensive  in  terms  of  cost  and  computation  time,  and  is 
thus  suitable  for  the  sub-model  that  is  integrated  into  3D 
simulation  code. 

(b)  The  ESP  model  provides  profile  information  on  concentra¬ 
tions  and  potentials  across  the  thickness  of  the  electrode. 
This  enables  heat  generation  rate  estimations  that  are  as 
accurate  as  the  Newman  model. 


Fig.  9.  Experimental  validation  of  the  ESP  model.  This  figure  compares  simulated 
discharged  curves  of  the  18650-type  lithium-ion  cell  (in-house  manufactured)  at 
various  rates.  The  lines  denote  model  predictions  while  circles  denote  the  experi¬ 
mental  data. 


These  two  advantages  are  essential  for  achieving  a  multi¬ 
dimensional  electrochemical-thermal  coupled  simulation 
required  for  lithium-ion  batteries. 
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Fig.  11.  A  schematic  of  the  two-way  electrochemical-thermal  coupled  simulation  method  of  lithium-ion  batteries. 


4.  Two-way  coupling  analysis 

4  A.  Numerical  procedure 

In  our  approach,  the  electrochemical  and  thermal  behaviors  of 
lithium-ion  batteries  are  simulated  using  two  separate  in-house 
codes.  Fig.  11  shows  a  schematic  of  our  two-way  electrochemical- 
thermal  coupled  simulation  method  for  lithium-ion  batteries. 

The  quasi-3D  porous  electrode  solver  is  applied  for  the  2D 
unrolled  plane  of  the  spirally  wound  electrodes.  Usually,  since  the 
curvature  radius  of  a  spirally  wound  object  is  sufficiently  large  (in 
relation  to  the  thickness  of  the  electrode  of  a  positive/negative  pole) 
even  if  it  treats  the  electrochemical  reaction  between  electrodes  at 
the  plane  locally,  the  potential  for  error  is  small  and  2D  deployment 
approximation  is  possible.  Thus,  the  ESP  model  could  be  imple¬ 
mented  into  each  2D  mesh.  Conservation  equations  of  charge  and 
lithium-ion  are  solved  in  the  2D  domain  by  estimating  their 
transport  across  the  electrode  thickness  using  the  ESP  model.  In  this 
meaning,  we  call  this  code  the  quasi-3D  solver.  Then,  potential  and 
lithium-ion  distributions  across  the  2D  unrolled  electrode  plane  are 
produced.  From  these  results,  we  can  obtain  the  distribution  of  the 
heat  generation  rate,  which  is  estimated  by  Eq.  (60). 

Next,  a  virtual  spirally  wound  electrode  is  constructed  by  the 
coordinate  transform  of  the  unrolled  electrode  plane  with  keeping 
the  heat  generation  rate  distribution  information.  Following  this, 
estimated  heat  generation  rate  data  are  mapped  to  the  3D  domain 
of  the  real  electrode  region,  which  is  used  in  the  3D  thermal  solver, 
and  the  3D  distribution  of  the  mapped  heat  generation  rate  qmapped 
is  reconstructed.  This  mapping  procedure  is  carried  out  by  means  of 
the  code  coupling  interface  (CCI)  [25  ,  which  is  supplied  by  AVL  List 
GmbH  and  which  makes  it  easy  for  both  codes  to  exchange  the 
specific  data. 


Overviews  of  the  3D  thermal  solver  and  the  coordinate  trans¬ 
form  in  the  data  mapping  procedure  are  presented  in  Appendix  A. 

The  3D  thermal  solver  simulates  the  temperature  distribution  of 
the  real  cell  using  the  mapped  heat  release  rate  arising  from  the 
electrochemical  reaction.  Next,  the  3D  mesh  of  the  spirally  wound 
electrode  region  in  the  real  cell  geometry  is  virtually  unrolled  onto 
the  2D  electrode  plane  by  the  reverse  coordinate  transform. 
Following  this,  temperature  data  are  mapped  to  the  2D  unrolled 
electrode  plane  by  means  of  CCI. 

The  abovementioned  two-way  data  exchange  process  between 
two  separate  solvers  is  performed  at  every  computational  time- 
step;  as  a  result,  both  electrochemical  and  thermal  behaviors  can 
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Fig.  12.  Comparison  between  time  evolutions  of  temperatures  for  the  18650-type 
lithium-ion  cell  obtained  from  experiment  and  those  predicted  by  the  two-way 
coupled  simulation  method.  The  time  evolutions  of  cell  potential  during  the  charge- 
discharge  cycle  are  also  shown  in  this  figure.  Applied  current  is  2.5  A  (5C),  and  ambient 
temperature  is  0  °C. 
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be  reproduced  simultaneously.  Additionally,  physical  properties  are 
updated  depending  on  the  temperature  at  every  computational 
time-step. 

4.2.  Implementation  of  ESP  model 

In  this  section,  the  formulation  necessary  to  implement  the  ESP 
model  into  the  quasi-3D  porous  electrode  solver  is  described  in 
detail. 

As  mentioned  above,  the  ESP  model  is  implemented  into  each 
2D  mesh  of  the  unrolled  electrode  plane.  However,  a  set  of  spirally 
wound  electrodes  consists  of  current  collectors,  composite  negative 
and  positive  electrodes,  and  a  separator.  Hence,  an  actual  unrolled 
electrode  plane  is  usually  in  the  order  of  100  pm  thick.  Accordingly, 
the  electrochemical  analysis  of  unrolled  electrodes  should  be 
strictly  performed  in  3D-mesh  with  an  extremely  large  aspect  ratio 
L/D  (L:  electrode  length,  D:  electrode  thickness). 

For  example,  the  3D  governing  equation  of  charge  conservation 
in  the  negative  electrolyte  phase  is  given  by 


Fig.  14.  Comparison  between  time  evolutions  of  temperatures  for  the  18650-type 
lithium-ion  cell  from  experiment  and  those  predicted  by  the  two-way  coupled  simu¬ 
lation  method.  The  time  evolutions  of  cell  potential  during  the  charge-discharge  cycle 
are  also  shown  in  this  figure.  Applied  current  is  2.5  A  (5C),  and  ambient  temperature  is 
20  °C. 
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where  a  x—y  plane  is  the  unrolled  electrode  plane.  In  the  quasi-3D 
porous  electrode  solver,  Eq.  (63)  reduces  to  the  2D  equation  with 
the  additional  source  term  as  follows: 
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The  last  term  on  the  left-hand  side  represents  the  charge 
transfer  from  the  anode  to  the  cathode  via  the  separator,  and  can  be 
formulated  using  the  ESP  model  as  follows: 
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Other  3D  governing  equations  for  charge  and  lithium-ion 
transfer  can  be  formulated  in  an  analogous  manner. 


5.  Temperature  measurement 

In  order  to  evaluate  the  calculation  accuracy  of  the  two-way 
electrochemical-thermal  coupled  simulation  method  presented  in 
this  paper,  temperatures  inside  the  cell  had  to  be  measured.  For  this 
purpose,  we  used  an  in-house  manufactured  18650-type  lithium- 
ion  cell  with  a  capacity  of  0.5  Ah.  The  cell  consists  of  Ni-based 
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Fig.  13.  Simulated  temperature  distribution  inside  the  cell  at  the  end  of  charge-  Fig.  15.  Simulated  temperature  distribution  inside  the  cell  at  the  end  of  charge- 

discharge  cycle.  Applied  current  is  2.5  A  (5C),  and  ambient  temperature  is  0  °C.  discharge  cycle.  Applied  current  is  2.5  A  (5C),  and  ambient  temperature  is  20  °C. 
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Fig.  16.  Comparison  between  time  evolutions  of  temperatures  for  the  18650-type 
lithium-ion  cell  from  experiment  and  those  predicted  by  the  two-way  coupled  simu¬ 
lation  method.  The  time  evolutions  of  cell  potential  during  the  charge-discharge  cycle 
are  also  shown  in  this  figure.  Applied  current  is  5  A  (IOC),  and  ambient  temperature  is 
20  °C. 


composite  oxides  positive  and  graphite  negative.  Experimental 
conditions  were  set  as  follows.  Ambient  temperatures  are  0  and 
20  °C.  The  charge-discharge  condition  is  Constant  Current  (CC) 
with  no  rest.  Applied  currents  are  2.5  and  5.0  A.  Cut-off  voltages  are 
4.1  V  for  charge  and  3.0  V  for  discharge. 

Two  18650-type  lithium-ion  cells  were  supplied  for  measure¬ 
ment  purposes,  while  type-T  (copper-constantan)  thermocouples 
were  used  for  the  measurements  themselves.  The  thermocouples 
were  0.13  mm  in  diameter.  Three  thermocouples  were  inserted  into 
one  cell  in  the  axial  direction,  and  another  three  were  inserted  into 
the  second  cell  along  the  radial  direction.  In  addition,  six  thermo¬ 
couples  were  placed  on  the  can  wall  of  each  cell. 

6.  Results  and  discussion 

6.1.  Comparisons  between  simulated  results  and  measured  results 

Fig.  12  shows  comparisons  between  the  time  evolutions  of 
temperatures  of  the  18650-type  lithium-ion  cell  obtained  via 
experiment  and  those  calculated  by  the  two-way  coupled 


time=2310  sec. 

Temp.  diff.  :  about  9‘ 

Fig.  17.  Simulated  temperature  distribution  inside  the  cell  at  the  end  of  charge- 
discharge  cycle.  Applied  current  is  5  A  (10C),  and  ambient  temperature  is  20  °C. 


simulation  method.  The  lower  two  sawtooth  lines  compare  time 
variations  of  cell  potential.  Applied  current  is  2.5  A  (5C).  The 
ambient  temperature  and  the  initial  temperature  of  the  prepared 
cell  were  set  at  0  °C.  It  should  be  noted  that  the  temperatures  in  this 
figure  reflect  the  average  value  of  the  measured  temperatures 
because  differences  at  each  measured  point  were  within  one  de¬ 
gree,  as  shown  in  Fig.  13.  However,  the  maximum  temperature 
difference  across  the  entire  cell  was  about  five  degrees  because  of 
the  large  heat  mass  at  the  upper  side.  It  is  noteworthy  that  the  time 
variations  of  temperature  measurements,  both  inside  the  cell  and 
on  the  can  wall,  were  simulated  within  three  degrees  of  accuracy, 
and  that  the  simulation  of  the  temperature  fluctuation  during  the 
charge-discharge  cycle  could  be  performed  with  high  accuracy. 

Other  validation  calculation  results  are  exhibited  in  the 
following  figures.  Fig.  14  shows  the  results  in  a  case  where  the 
applied  current  is  2.5  A  (5C)  and  ambient  temperature  is  20  °C. 
From  Fig.  15,  it  can  be  seen  that  the  maximum  temperature  dif¬ 
ference  across  the  entire  cell  becomes  about  four  degrees  at  the  end 
of  the  charge-discharge  cycle  in  this  case.  Additionally,  Fig.  16 
shows  the  result  of  a  case  where  the  applied  current  is  5  A  (10C) 
and  ambient  temperature  is  20  °C,  in  which  case  the  maximum 
temperature  difference  across  the  entire  cell  reaches  about  nine 
degrees,  as  shown  in  Fig.  17. 

Based  on  these  validation  results,  it  can  be  seen  that  the  two- 
way  electrochemical-thermal  coupled  simulation  method  pre¬ 
sented  in  this  paper  is  capable  of  simulating  temperatures  inside 
the  cell  and  on  the  can  wall  with  an  accuracy  level  within  a  few 
degrees. 


6.2.  Contribution  of  each  heat  release  source  term 

In  this  section,  we  will  devote  a  little  more  space  to  discussing 
the  contribution  of  each  heat  source  term  to  total  electrochemical 
heat  release  rate. 

As  previously  mentioned,  the  local  electrochemical  heat  release 
rate  consists  of  three  heat  source  terms  that  are  expressed  in  Eq.  (1 ) 
and/or  Eq.  (60): 


Qelchem  —  Qohmic  "I"  Qirrev  "T  Qrev  (66) 

The  first  term  on  the  right-hand  side  refers  to  the  ohmic  heat 
caused  by  ionic  resistance.  The  second  term  represents  the  irre¬ 
versible  heat  arising  from  overpotential  and  film  resistance.  The 
third  term  expresses  the  reversible  heat  resulting  from  the  entropic 
effect  of  lithium  intercalation  or  deintercalation.  These  three  terms 
are  formulated  as  follows: 

Qohmic  =  j)  +  5^  (^f^^e,  j '^^e,  j 

j  =  n,  p  j  =  n,  sep,  p 

+  kD ,  ce ,  j '  j) 

(67) 

Qirrev  =  ^  ]  ^s,  j*n,  j  (^s,  j  ~  fie,  j  ~  (68) 

j  =  n,  p 

-  dU; 

Qrev  =  as  jln  J-T—  (69) 

j  =  n,  p 

Fig.  18  displays  the  effect  of  each  three  heat  source  term  on  time 
variations  of  temperature  of  the  cell  during  the  charge-discharge 
cycle  at  the  5C  rate  (2.5  A),  0  °C.  Here,  only  one-cycle  discharge- 
charge  process  is  shown.  As  indicated  in  Fig.  18(a),  a  small  tem¬ 
perature  fluctuation  is  observed  inside  the  cell.  Fig.  18(b)  and  (c) 
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Fig.  18.  Effect  of  each  heat  source  term  on  time  variations  of  temperature  of  the  cell  during  the  charge-discharge  cycle  at  the  5C  rate  (2.5  A),  0  °C.  The  small  fluctuation  of  the 
temperature  behavior  is  mainly  caused  by  the  reversible  heat  resulting  from  entropic  effect  due  to  lithium-ion  intercalation  or  deintercalation. 
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(b)  Distributed  pattern  of  the  total  electrochemical  heat  release  rate 
at  DOD=30%. 
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(c)  Distributed  pattern  of  the  total  electrochemical  heat  release  rate 
at  DOD=50%. 


(a)  Charge  and  discharge  curves. 


Fig.  19.  Time  evolutions  of  distributed  patterns  of  the  total  electrochemical  heat  release  rate  across  the  entire  electrode  plane  during  discharge  process  at  the  10C  rate  (5  A),  20  °C. 


shows  the  time  variations  of  the  contribution  ratio  for  each  of  the 
three  heat  source  terms  in  relation  to  the  total  electrochemical  heat 
release  rate  during  the  discharge  and  charge  process,  respectively. 
From  these  figures,  it  is  clear  that  the  ohmic  heat  is  similar  in  value 
between  discharge  and  charge  and  remains  constant  over  the 
entire  process.  In  addition,  it  can  be  seen  that  the  irreversible  heat 
is  also  similar  in  value  between  both  processes  and  maintains 
nearly  constant.  As  a  result,  a  remarkable  difference  is  observed  in 
the  reversible  heat,  which  is  expressed  by  Eq.  (69). 

The  time  variation  of  the  total  electrochemical  heat  release  rate 
is  directly  reflected  in  the  thermal  behavior  of  the  cell,  as  shown  in 
Fig.  18.  Furthermore,  the  reversible  heat  governs  the  time  variation 
behavior  of  the  heat  generation  rate.  Accordingly,  these  results  lead 
to  the  conclusion  that  the  small  fluctuation  noted  in  the  tempera¬ 
ture  behavior  is  caused  primarily  by  the  reversible  heat. 

Reversible  heat  variations  depend  on  the  local  state  of  charge 
(SOC)  of  the  anode  and  cathode  active  materials  during  the 
discharge  and  charge  process.  In  order  to  evaluate  the  reversible 
heat,  the  entropy  change  AS  value  of  the  active  materials,  which  is 
related  to  the  derivative  of  the  open  circuit  potential  (OCP)  with 
respect  to  temperature,  is  necessary.  In  this  study,  instead  of 
measuring  AS  directly,  we  used  a  model  formulation  that  predicts 


the  entropy  change  based  on  experimental  data  of  the  OCP.  Details 
of  the  model  used  for  entropy  change  are  provided  in  Appendix  B. 

6.3.  Distribution  of  heat  generation  rate 

In  the  two-way  coupled  simulation  method  presented  in  this 
paper,  the  electrochemical  phenomena  of  lithium-ion  batteries  are 
simulated  using  the  quasi-3D  porous  electrode  solver.  Hence,  the 
2D  distributed  information  of  a  heat  generation  rate  arising  from 
electrochemical  reactions  across  an  entire  electrode  plane  can  be 
obtained. 

Fig.  19  shows  the  distribution  of  a  heat  generation  rate  across  an 
entire  unrolled  electrode  plane.  Fig.  19(a)  shows  comparisons  be¬ 
tween  charge  and  discharge  curves  obtained  from  experiments  and 
simulated  results  at  the  10C  rate  (5  A),  20  °C,  while  Fig.  19(b)  and  (c) 
shows  distributed  patterns  of  the  total  electrochemical  heat  gen¬ 
eration  rate  when  the  DOD  in  the  discharge  process  equals  30%  and 
50%,  respectively. 

At  DOD  =  30%,  the  maximum  heat  generation  rate  puts  in  an 
appearance  at  the  end  near  the  negative  tab,  while  the  minimum 
value  occurs  at  the  central  region.  In  contrast,  the  distributed 
pattern  of  the  heat  generation  rate  turns  around  at  DOD  =  50%.  In 
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other  words,  the  distributed  pattern  of  the  total  electrochemical 
heat  generation  rate  varies  moment  by  moment  depending  on  the 
DOD  or  SOC. 

From  these  figures,  it  can  be  seen  that  the  range  between 
maxima  and  minima  of  the  total  electrochemical  heat  generation 
rate  across  the  entire  electrode  plane  reaches  5%  at  the  IOC  rate 
discharge  in  the  case  of  18650-type  cell  of  0.5  Ah.  As  for  the  5C  rate 
discharge,  it  reaches  the  vicinity  of  1%;  however,  the  figure  showing 
the  calculation  results  is  omitted  here.  Accordingly,  it  is  easy  to 
imagine  that  the  ununiformity  of  the  total  heat  generation  rate 
across  the  entire  electrode  plane  becomes  larger  as  the  applied 
current  rate  increases. 

7.  Conclusions 

An  enhanced  single  particle  (ESP)  charge-discharge  model  for 
lithium-ion  batteries  has  been  developed.  This  ESP  model  handles 
each  negative  and  positive  electrode  as  a  single  spherical  particle  in 
the  electrolyte  phase.  The  potential  and  lithium-ion  concentration 
in  the  electrolyte  phase  are  set  at  the  representative  position  in 
each  electrode,  and  the  distributions  of  these  physical  properties 
across  the  thickness  of  the  electrode  are  approximated  by  parabolic 
profiles.  The  use  of  these  model  features  facilitate  more  accurate 
calculations  of  the  charge-discharge  curve  under  high  rate  condi¬ 
tions  than  are  possible  using  the  conventional  SP  model  because 
the  potential  drop  in  the  electrolyte  phase  can  be  simulated  in  the 
new  model. 

This  ESP  model  is  inexpensive  in  terms  of  cost  and  computation 
time  requirements  because  it  uses  the  mass-point  system.  In 
addition,  because  the  accurate  estimation  of  the  heat  generation 
rate  requires  information  on  potential  and  lithium-ion  concentra¬ 
tion  in  the  solution  phase,  the  ESP  model  is  capable  of  estimating 
the  heat  generation  rate  caused  by  electrochemical  reactions  more 
accurately  than  previous  methods. 

The  above  two  advantages  have  contributed  significantly  to  the 
development  of  our  multi-dimensional  two-way  electrochemical- 
thermal  coupled  simulation  method.  This  method  has  been 
applied  to  a  thermal  behavior  analysis  of  the  charge-discharge 
cycle  of  an  18650-type  lithium-ion  cell  battery,  where  good 
agreement  between  the  simulation  and  actual  measurements  was 
achieved.  Furthermore,  the  results  clearly  show  that  the  small 
temperature  behavior  fluctuation  observed  during  the  charge- 
discharge  cycle  is  caused  primarily  by  the  reversible  heat  resulting 
from  the  entropic  effect  arising  from  lithium  intercalation  or 
deintercalation. 


wound  electrode  is  in  a  state  where  the  anode,  the  cathode,  and 
separator  are  impregnated  with  the  electrolyte,  the  space  between 
the  positive  and  negative  electrodes  is  narrow  and  temperature 
differences  are  small.  Furthermore,  since  it  is  difficult  to  treat  them 
individually  when  considering  the  spatial  resolution  of  3D  thermal 
analysis,  an  electrode  is  treated  as  a  single  unit.  At  that  time,  the 
average  physical  properties  value  was  given  as  follows. 

First,  the  individual  thermo  physical  properties  values  of  the 
electrolyte,  anode,  cathode,  and  separator  are  specified.  In  the  pores 
of  the  positive/negative  electrodes  and  separator,  which  are  filled 
with  electrolyte,  the  value  of  each  physical  property  is  estimated  as 
a  volume  average  of  the  solid  and  liquid  phases  using  their  per¬ 
centage  of  the  void.  Finally,  the  average  physical  properties  value  of 
a  laminated  electrode  is  computed.  The  density  and  specific  heat 
are  estimated  as  a  weight-average  of  the  thickness  of  each  part  for 
electrode  laminated  structure.  On  the  other  hand,  while  the  ther¬ 
mal  conductivity  of  a  radial  (electrode  thickness)  direction  differs 
from  those  of  the  other  two  directions  (a  circumference,  an  axis), 
they  are  given  by  the  following  formulas  for  all  directions. 


E  ¥t 

At’z  =  W  Ar 


(A2) 


Here,  the  subscripts  t,  z  and  r  refer  to  the  circumferential,  axis 
and  radial  directions  in  the  cylindrical  coordinates  system  based  on 
the  central  axis  of  a  spirally  wound  object.  Additionally,  <5  is  thick¬ 
ness,  and  the  subscript  i  shows  the  lamination  component  of 
electrode  for  one  cycle. 

Next,  the  coordinate  transformation  for  data  mapping  between 
a  2D  electrode  plane  mesh  and  a  3D  spirally  wound  thermal  model 
mesh  is  described.  Here,  the  relationship  between  the  x-coordinate 
of  the  2D  electrode  deployment  direction  and  the  radius  r  and  an 
angle  6  of  a  spirally  wound  object  are  as  follows: 


r 


r0  + 


(A3) 


0 


2irr0 

a 


(A4) 


Here,  vq  is  the  inner  end  radius  and  a  is  a  pitch  of  one  cycle  of 
electrode  lamination  as  shown  in  Fig.  20.  It  should  be  noted  that  0 
[rad]  is  not  2tt  cycle  but  monotonically  increases  with  the  number 
of  windings. 


Appendix  A.  Cell  thermal  solver  and  coordinate 
transformation  for  data  mapping  between  electrode  and  cell 
thermal  models 

In  the  thermal  analysis  of  a  battery  cell  interior,  supposing  the 
convective  velocities  of  gas  enclosed  in  a  cell  and  an  electrolyte  are 
minute  and  can  be  ignored,  the  governing  equation  of  a  thermal 
field  turns  into  the  following  3D  heat  conduction  equation. 

pCpg  =  V(AV7)  +  q  (Al) 

Here,  p,  Cp  and  A  are  density,  specific  heat  and  thermal  con¬ 
ductivity,  respectively. 

A  normal  battery  cell  consists  of  a  top  cap,  packing,  a  collector 
tab,  an  insulator,  a  spacer  (gas),  etc.,  in  addition  to  the  electrode  and 
the  can.  The  thermo  physical  properties  value  (density,  specific 
heat,  thermal  conductivity)  of  each  of  those  constructional  ele¬ 
ments  is  provided  for  in  the  simulation.  However,  although  a 


Appendix  B.  Model  for  entropy  change  during  intercalation 
and  deintercalation 

In  this  Appendix,  we  describe  the  model  used  to  predict  the 
entropy  change  necessary  for  estimating  the  reversible  heat  (see 
Section  6.2).  The  entropy  change  due  to  intercalation  and  dein¬ 
tercalation  is  experimentally  obtained  by  measuring  the 


Fig.  20.  Illustration  of  a  spirally  wound  object. 
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Fig.  21.  Comparison  of  the  model  prediction  and  experimental  results,  (a)  Open  circuit  potential  U.  (b)  Entropy  change  (5 U/dT)  =  kS/F.  The  solid  line  indicates  the  present  results,  the 
dashed  line  indicates  the  model  prediction  by  Wong  and  Newman  [29],  and  the  symbol  indicates  the  experimental  results  provided  by  Thomas  et  al.  [28]. 


temperature  dependence  of  the  open  circuit  potential  (OCP). 
However,  performing  this  measurement  is  difficult  because  it  takes 
a  significant  amount  of  time  to  reach  the  equilibrium  state  (at 
which  the  OCP  can  be  measured)  after  the  temperature  is  changed. 
In  contrast,  the  model  described  in  this  Appendix  predicts  the 
temperature  dependence  based  solely  on  the  OCP  data  at  a  fixed 
temperature. 

The  model  is  premised  on  the  assumption  that  the  active  ma¬ 
terial  consists  of  N  independent  components  that  have  different 
reaction  properties.  This  is  a  model  for  describing  the  stepwise 
change  of  the  OCP  depending  on  the  SOC,  which  is  used  by  Zhang 
et  al.  [26]  to  derive  a  kinetic  model. 

We  begin  by  introducing  several  variables  into  the  model 
description.  Let  0  be  the  ratio  of  the  number  of  intercalated  lithium 
atoms  to  the  total  number  of  the  sites  in  the  active  material.  In 
addition,  let  Ok  and  Ok  denote  the  ratios  of  the  number  of  the  sites, 
and  the  number  of  lithium  atoms  in  kth  component,  respectively,  to 
the  total  number  of  sites  in  the  active  material: 

N  N 

E0'<  =  1’  =  (B1) 

k  =  1  k=  1 

Following  the  discussion  in  Ref.  26],  starting  from  the  Butler- 
Volmer  kinetics  for  multi-component  material,  the  following  rela¬ 
tion,  which  holds  at  the  equilibrium  state,  is  derived: 


U 


U°k  + 


(B2) 


where  U®  is  a  quantity  related  to  the  reaction  rate,  and  ?/<  is  the 
order  of  the  reaction  of  the  kth  component.  Solving  Eq.  B2)  with 
respect  to  Ok  and  substituting  into  Eq.  (Bl)  yields 
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This  relationship  between  0  and  U  with  the  model  parameters 
Ok,  U®,  and  ?/<  coincides  with  the  OCP  model  proposed  semi- 
empirically  by  Verbrugge  and  Koch  in  their  work  on  modeling 
graphite  electrodes  [27].  Differentiation  of  Eq.  (B3)  with  respect  to  T 
gives  the  quantity  related  to  the  entropy  change  AS  through  (dU/ 
0T)  =  A  S/F: 


where 


£‘"'xp(pr(u-<)' 


(B4) 


(B5) 


Once  we  have  the  (B3)  relationship  fitted  to  a  0- U  curve 
measured  experimentally,  the  predicted  values  of  entropy  change 
are  drawn  immediately  using  the  same  set  of  model  parameters  Ok, 
U°,  and  ffe.  Note  that  the  coefficient  which  is  related  to  the  re¬ 
action  rate,  is  intrinsically  a  function  of  temperature.  In  the  above 
derivation,  however,  we  ignore  the  dependency  of  on  temper¬ 
ature  in  view  of  the  fact  that  the  dependency  vanishes  in  cases 
where  the  frequency  factor  value  is  common  to  the  anodic  and 
cathodic  reaction. 

The  adequacy  of  the  model  is  assessed  by  comparison  with  the 
experimental  results  by  Thomas  et  al.  [28]  They  measured  the  OCP 
and  dU/dT  of  lithium  manganese  spinel  (LiyMn204)  at  various 
values  of  0.  The  values  of  the  model  parameters  Ok,  U®,  and  ?/<  were 
determined  by  fitting  Eq.  (B3)  with  N  =  6  to  the  experimental  data 
of  the  OCP.  The  method  we  employed  to  accomplish  this  is  the 
standard  nonlinear  general  reduced  gradient  (GRG)  optimization 
algorithm.  We  then  estimated  the  values  of  dU/dT  using  Eq.  (B4). 
Fig.  21  compares  the  results  of  the  present  model  and  Ref.  27. 
Although  the  present  model  underestimates  the  value  of  (dU/ 
9T)  =  A  S/F  in  the  intermediate  region  of  0 ,  the  overall  entropy 
change  trend  is  captured  well.  In  the  figure,  the  results  of  the 
model  provided  by  Wong  and  Newman  [29]  are  also  shown.  In 
their  model,  both  the  OCP  and  dU/dT  are  predicted.  The  accuracy  of 
the  entropy  change  is  almost  the  same  as  the  present  model, 
whereas  the  reproduction  of  the  OCP  of  the  present  model  is  more 
accurate. 


List  of  symbols 

as  specific  interfacial  area  of  an  electrode  (m2  nrT3) 

ce  Li+  concentration  in  the  electrolyte  phase  (mol  m  3) 

cs  Li  concentration  in  the  solid  particles  (mol  m-3) 

cse  Li  concentration  at  the  surface  of  the  solid  particles 

(mol  m-3) 
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De  diffusion  coefficient  of  Li+  in  the  electrolyte  phase 

(m2  s_1) 

Ds  diffusion  coefficient  of  Li+  in  the  solid  phase  (m2  s_1) 

F  Faraday’s  constant  (96.487C  mol-1) 

7ap p  applied  current  density  (A  rrT2) 

i0  exchange  current  density  (A  nrT2) 

in  superficial  current  density  (A  rrT2) 

j  local  volumetric  transfer  current  density  due  to  charge 

transfer  (A  rrT3) 

lse  diffusion  length  of  Li+  from  solid-electrolyte  interface 

into  solid  phase  (m) 

L  thickness  of  n,  sep  or  p  (m) 

q  heat  generation  rate  (W  nrT3) 

R  universal  gas  constant  (8.3143  J/(mol  1  •  I<)) 

Rf  film  resistance  on  an  electrode  surface  (Q  m2) 

r  radial  coordinate  (m) 

rs  radius  of  the  spherical  particles  (m) 

T  absolute  temperature  (I<) 

t  time  (s) 

t°+  transference  number  of  Li+  in  solution 

U  open-circuit  potential  of  an  electrode  reaction  (V) 

x,  y  coordinates  on  the  unrolled  electrode  plane  (m) 
z  coordinate  across  the  thickness  of  the  electrode  (m) 

Greek  symbols 

aa,ac  anodic  and  cathodic  transfer  coefficients  for  an  electrode 
reaction 

£  volume  fraction  of  a  phase 

p  surface  overpotential  of  an  electrode  reaction  (V) 

k  conductivity  of  the  electrolyte  (S  m-1) 

kd  diffusional  conductivity  of  the  electrolyte  (Am-1) 

a  conductivity  of  the  electrode  (S  itT1) 

(ps  electrical  potential  in  the  solid  phase  (V) 

4>e  electrical  potential  in  the  electrolyte  phase  (V) 

Subscript  and  superscript 
eff  effective 


j  n  or  p 

Li  lithium  species 

n  negative  electrode 

p  positive  electrode 

sep  separator 
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